{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "SzKwuqYESWwm"
   },
   "source": [
    "##### Copyright 2022 The Cirq Developers"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "cellView": "form",
    "id": "4yPUsdJxSXFq"
   },
   "outputs": [],
   "source": [
    "# @title Licensed under the Apache License, Version 2.0 (the \"License\");\n",
    "# you may not use this file except in compliance with the License.\n",
    "# You may obtain a copy of the License at\n",
    "#\n",
    "# https://www.apache.org/licenses/LICENSE-2.0\n",
    "#\n",
    "# Unless required by applicable law or agreed to in writing, software\n",
    "# distributed under the License is distributed on an \"AS IS\" BASIS,\n",
    "# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.\n",
    "# See the License for the specific language governing permissions and\n",
    "# limitations under the License."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "b5cbe21a91f4"
   },
   "source": [
    "# Coherent vs Incoherent Error with XEB"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "zC1qlUJoSXhm"
   },
   "source": [
    "<table class=\"tfo-notebook-buttons\" align=\"left\">\n",
    "  <td>\n",
    "    <a target=\"_blank\" href=\"https://quantumai.google/cirq/noise/qcvv/coherent_vs_incoherent_xeb>\"><img src=\"https://quantumai.google/site-assets/images/buttons/quantumai_logo_1x.png\" />View on QuantumAI</a>\n",
    "  </td>\n",
    "  <td>\n",
    "    <a target=\"_blank\" href=\"https://colab.research.google.com/github/quantumlib/Cirq/blob/main/docs/noise/qcvv/coherent_vs_incoherent_xeb.ipynb\"><img src=\"https://quantumai.google/site-assets/images/buttons/colab_logo_1x.png\" />Run in Google Colab</a>\n",
    "  </td>\n",
    "  <td>\n",
    "    <a target=\"_blank\" href=\"https://github.com/quantumlib/Cirq/blob/main/docs/noise/qcvv/coherent_vs_incoherent_xeb.ipynb\"><img src=\"https://quantumai.google/site-assets/images/buttons/github_logo_1x.png\" />View source on GitHub</a>\n",
    "  </td>\n",
    "  <td>\n",
    "    <a href=\"https://storage.googleapis.com/tensorflow_docs/Cirq/docs/noise/qcvv/coherent_vs_incoherent_xeb.ipynb\"><img src=\"https://quantumai.google/site-assets/images/buttons/download_icon_1x.png\" />Download notebook</a>\n",
    "  </td>\n",
    "</table>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "dd95be2a71eb"
   },
   "source": [
    "This notebook demonstrates how to use Cross-Entropy Benchmarking (XEB) end-to-end to compare coherent gate parameter error and incoherent depolarization error. It will mimic a small device graph of `cirq.GridQubit` pairs and simulate two-qubit XEB benchmarking circuits on them with noise models that introduce coherent and incoherent error before visualizing and comparing the results.\n",
    "\n",
    "For more information on types of error, see [Average, Pauli and Incoherent Error](../../google/calibration.md#average-pauli-and-incoherent_error). \n",
    "\n",
    "For more information on noise and noise models, see [Noise](../).\n",
    "\n",
    "For more information on how XEB works in detail, see [XEB Theory](./xeb_theory.ipynb) notebook."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "57flCRpJpup3"
   },
   "source": [
    "## Setup"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "bd9529db1c0b"
   },
   "outputs": [],
   "source": [
    "try:\n",
    "    import cirq\n",
    "except ImportError:\n",
    "    print(\"installing cirq...\")\n",
    "    !pip install --quiet cirq\n",
    "    print(\"installed cirq.\")\n",
    "    import cirq\n",
    "\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "ace31cc4d258"
   },
   "source": [
    "## Generate Random Circuits\n",
    "\n",
    "First create a circuit library of 20 random, two-qubit circuits which uses `cirq.SQRT_ISWAP` as the entangling two-qubit gate, using `cirq.experiments.random_quantum_circuit_generation.generate_library_of_2q_circuits\n",
    "`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "6fa8a0266b7a"
   },
   "outputs": [],
   "source": [
    "from cirq.experiments import random_quantum_circuit_generation as rqcg\n",
    "\n",
    "RANDOM_SEED = np.random.RandomState(53)\n",
    "\n",
    "circuit_library = rqcg.generate_library_of_2q_circuits(\n",
    "    n_library_circuits=20, two_qubit_gate=cirq.SQRT_ISWAP, random_state=RANDOM_SEED\n",
    ")\n",
    "max_depth = 100\n",
    "cycle_depths = np.arange(3, max_depth, 20)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "2e0f9de60ef1"
   },
   "source": [
    "## Prepare Noise Models\n",
    "\n",
    "To compare coherent and incoherent error, use two different noise models, which introduce coherent and incoherent error respectively. Also included are a noiseless model and a \"fused\" model that combines both coherent and incoherent error, for comparison. See [Noise](../) for the difference between coherent and incoherent error. \n",
    "\n",
    "The coherent noise model is represented by a perturbed `SQRT_ISWAP` gate, which implements a slightly different unitary than a true `SQRT_ISWAP`. Perturbed gates are caused by hardware losing its calibrated tuning, which typically occurres naturally over time between device calibrations. \n",
    "\n",
    "The incoherent noise model to be compared with is a `cirq.DepolarizingChannel`, representing a low chance (`5e-3`) for the state of each qubit to depolarize after an operation is performed on it, as if a random Pauli operator was applied to it.\n",
    "\n",
    "An noiseless sampler simply has no noise. \n",
    "\n",
    "Finally, the `fused_sampler` adds both coherent and incoherent error to each moment of the circuit, applied in order, according to the `noisy_moment` function of the `FusedNoiseModel` class. It uses the same noise models as the coherent-only and incoherent-only models, perturbation and depolarization, but applies both in sequence. \n",
    "\n",
    "Start by creating the necessary elements to create the coherent and fused samplers. The following code shows how `SQRT_ISWAP` can be written as a specific `cirq.PhasedFSimGate`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "a598f743d18a"
   },
   "outputs": [],
   "source": [
    "# build sqrt_iswap gate from a phased fsim gate\n",
    "sqrt_iswap_as_phased_fsim = cirq.PhasedFSimGate.from_fsim_rz(\n",
    "    theta=-np.pi / 4, phi=0, rz_angles_before=(0, 0), rz_angles_after=(0, 0)\n",
    ")\n",
    "\n",
    "# check that the unitaries are the same\n",
    "np.testing.assert_allclose(\n",
    "    cirq.unitary(sqrt_iswap_as_phased_fsim), cirq.unitary(cirq.SQRT_ISWAP), atol=1e-8\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "31492475ce1b"
   },
   "source": [
    "Next, create a perturbed version of `SQRT_ISWAP`. Note the phi angle uniformly chosen in the range $[0, 2/16*\\pi]$, which was previously `0`. Additionally, create the gate substitution function, `_sub_iswap`, which replaces `SQRT_ISWAP` gates in the simulated circuit with the perturbed version."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "da1de99252fa"
   },
   "outputs": [],
   "source": [
    "from collections import defaultdict\n",
    "\n",
    "phi_angles = defaultdict(lambda: RANDOM_SEED.uniform(low=0, high=2 / 16 * np.pi))\n",
    "\n",
    "\n",
    "def _sub_iswap(op):\n",
    "    if op.gate == cirq.SQRT_ISWAP:\n",
    "        # add coherent error to the phi angle while creating a sqrt iswap gate\n",
    "        perturbed_sqrt_iswap = cirq.PhasedFSimGate.from_fsim_rz(\n",
    "            theta=-np.pi / 4,\n",
    "            phi=phi_angles[op.qubits],\n",
    "            rz_angles_before=(0, 0),\n",
    "            rz_angles_after=(0, 0),\n",
    "        )\n",
    "\n",
    "        return perturbed_sqrt_iswap.on(*op.qubits)\n",
    "    return op"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "hdaQLuUyX78g"
   },
   "source": [
    "For the \"fused\" model, create a `FusedNoiseModel` class, which applies multiple noise models in sequence to each moment. This will create a noise model which incorporates both coherent and incoherent noise. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "OKsBPXehPYxH"
   },
   "outputs": [],
   "source": [
    "from typing import Sequence\n",
    "\n",
    "\n",
    "class FusedNoiseModel(cirq.NoiseModel):\n",
    "    \"\"\"A noise model that applies other noise models in sequence\"\"\"\n",
    "\n",
    "    def __init__(self, models: Sequence[cirq.NoiseModel]):\n",
    "        self.models = models\n",
    "\n",
    "    def noisy_moment(self, moment: 'cirq.Moment', system_qubits: Sequence['cirq.Qid']):\n",
    "        \"\"\"Produce a list of moments by applying each model in sequence to the available moments\"\"\"\n",
    "\n",
    "        moments = [moment]\n",
    "        for model in self.models:\n",
    "            new_moments = []\n",
    "            for moment in moments:\n",
    "                ret = model.noisy_moment(moment, system_qubits)\n",
    "                # handle constituent noise models returning lists of or singleton moments\n",
    "                new_moments.extend([ret] if isinstance(ret, cirq.Moment) else ret)\n",
    "            moments = new_moments\n",
    "        return moments"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "BjcqjO95X6ny"
   },
   "source": [
    "Finally, initialize the noise model objects: \n",
    "\n",
    "1. For coherent error, use the `_sub_iswap` function to create a `cirq.GateSubstitutionNoiseModel`, which replaces the `SQRT_ISWAP` gate with the perturbed version. \n",
    "\n",
    "2. For incoherent error, create a `cirq.ConstantQubitNoiseModel` which depolarizes every qubit in the same way. \n",
    "\n",
    "3. The noiseless sampler will have a noise of `None`.\n",
    "\n",
    "4. Build the fused noise model with the other two from the `FusedNoiseModel` class, applying each constituent noise model in sequence.\n",
    "\n",
    "... and create the samplers with `cirq.DensityMatrixSimulator` for each noise model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "kWdgExibq8JU"
   },
   "outputs": [],
   "source": [
    "# create noise models\n",
    "coherent_noise = cirq.devices.noise_model.GateSubstitutionNoiseModel(_sub_iswap)\n",
    "incoherent_noise = cirq.ConstantQubitNoiseModel(cirq.depolarize(5e-3))\n",
    "noiseless = None\n",
    "fused_noise = FusedNoiseModel([coherent_noise, incoherent_noise])\n",
    "\n",
    "# create samplers by passing noise models to simulator constructor\n",
    "coherent_sampler = cirq.DensityMatrixSimulator(noise=coherent_noise)\n",
    "incoherent_sampler = cirq.DensityMatrixSimulator(noise=incoherent_noise)\n",
    "noiseless_sampler = cirq.DensityMatrixSimulator(noise=noiseless)\n",
    "fused_sampler = cirq.DensityMatrixSimulator(noise=fused_noise)\n",
    "\n",
    "# from here on out, all experiment steps will be performed on\n",
    "#   samplers and their results in the following order.\n",
    "samplers = [coherent_sampler, incoherent_sampler, noiseless_sampler, fused_sampler]\n",
    "noise_titles = ['Coherent', 'Incoherent', 'Noiseless', 'Fused']"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6fbec4054504"
   },
   "source": [
    "## Mock a device topology\n",
    "\n",
    "In order to more closely mimic XEB on real hardware, create a graph of qubits to mimic a small device from a couple `cirq.GridQubit`s. XEB benchmark circuits will be created to test multiple qubit pairs simultaneously on this graph, without overlap. \n",
    "\n",
    "In this example, only five of the six created qubits are used in order to produce an irregular/non-rectangular graph, and simplify simulation by using fewer qubits."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "af1e12169e88"
   },
   "outputs": [],
   "source": [
    "import networkx as nx\n",
    "import itertools\n",
    "\n",
    "# a set of six example GridQubits\n",
    "qubits = cirq.GridQubit.rect(3, 2, 4, 3)\n",
    "# only use the first num_qubits many qubits\n",
    "num_qubits = 5\n",
    "qubits = qubits[:num_qubits]\n",
    "\n",
    "# create graph from adjacent qubits\n",
    "graph = nx.Graph((q1, q2) for (q1, q2) in itertools.combinations(qubits, 2) if q1.is_adjacent(q2))\n",
    "pos = {q: (q.row, q.col) for q in qubits}\n",
    "nx.draw_networkx(graph, pos=pos)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "5d6f40acc7bf"
   },
   "source": [
    "Set up the active qubit pair combinations for the sampler function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "fbaf9c5c5335"
   },
   "outputs": [],
   "source": [
    "combs_by_layer = rqcg.get_random_combinations_for_device(\n",
    "    n_library_circuits=len(circuit_library),\n",
    "    n_combinations=10,\n",
    "    device_graph=graph,\n",
    "    random_state=RANDOM_SEED,\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "423f56a942f6"
   },
   "source": [
    "## Fidelity Data Collection\n",
    "### Sample Circuits\n",
    "\n",
    "The following call will zip together the `circuit_library` circuits according to `combs_by_layer` into larger circuits that evaluate multiple qubit pairs in parallel, before sampling them. Remember that `samplers` contains coherent and incoherent noise samplers, an errorless sampler, and a fused sampler. Mapping `sample_2q_xeb_circuits` to the list of samplers produces sampled data, generated with identical parameters except for the difference in noise model. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "aab08193f046"
   },
   "outputs": [],
   "source": [
    "import cirq.experiments.xeb_sampling as xeb_sampling\n",
    "import cirq.experiments.xeb_fitting as xeb_fitting\n",
    "\n",
    "\n",
    "# use the same circuit library, cycle depths, combinations by layer,\n",
    "#   random seed and repetitions for both noise models' samplers.\n",
    "def sample_df(sampler):\n",
    "    return xeb_sampling.sample_2q_xeb_circuits(\n",
    "        sampler=sampler,\n",
    "        circuits=circuit_library,\n",
    "        cycle_depths=cycle_depths,\n",
    "        combinations_by_layer=combs_by_layer,\n",
    "        shuffle=RANDOM_SEED,\n",
    "        repetitions=1000,\n",
    "    )\n",
    "\n",
    "\n",
    "sampled_dfs = list(map(sample_df, samplers))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "16b5a4a86384"
   },
   "source": [
    "Compute the circuit fidelities from the sampled probabilities for each sampled dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "8dd3ed7024ef"
   },
   "outputs": [],
   "source": [
    "# use the same circuit library and cycle_depths for both\n",
    "#   noise models' sampled data.\n",
    "def estimate_fidelities(sampled_df):\n",
    "    return xeb_fitting.benchmark_2q_xeb_fidelities(\n",
    "        sampled_df=sampled_df, circuits=circuit_library, cycle_depths=cycle_depths\n",
    "    )\n",
    "\n",
    "\n",
    "circuit_fidelities = list(map(estimate_fidelities, sampled_dfs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QLxjg-ZPzZsJ"
   },
   "source": [
    "Estimate the by-layer fidelities from the circuit fidelities for circuits of varying lengths for each set of fidelities."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "05748707cacc"
   },
   "outputs": [],
   "source": [
    "fitted_fidelities = list(map(xeb_fitting.fit_exponential_decays, circuit_fidelities))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "QDPtNQlPz3BL"
   },
   "source": [
    "## Visualizing Fidelity by Noise Model\n",
    "\n",
    "### Qubit Pair Heatmaps\n",
    "The following cell visualizes the by-cycle fidelity error in qubit pair heatmaps."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "WJbgW4r1Guuf"
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "fig, axes = plt.subplots(2, 2, figsize=(13, 13))\n",
    "for fidelities, ax, title in zip(fitted_fidelities, axes.flat, noise_titles):\n",
    "    # pull out the by-layer/cycle fidelity and turn it into error with 1-fidelity.\n",
    "    heatmap_data = {\n",
    "        pair: (1.0 - layer_fid) for ((_, _, pair), layer_fid) in fidelities.layer_fid.items()\n",
    "    }\n",
    "    cirq.TwoQubitInteractionHeatmap(heatmap_data, vmin=0, vmax=0.025).plot(ax)\n",
    "\n",
    "    ax.set_title(title)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "B-wy6SdblptC"
   },
   "source": [
    "The heatmaps reveal that the estimated error is quite similar when compared across qubits in the mock device. \n",
    "\n",
    "The coherent error heatmap shows the most variation in qubit pair fidelity error, though all values are within `0.003` of each other. The perturbed `SQRT_ISWAP` gate produces a consistent change in angle that interacts differently with the random single-qubit rotations added in the two-qubit circuits of the `circuit_library`, causing this variance.\n",
    "\n",
    "The incoherent error heatmap shows little to no variation in error. This makes sense because incoherent depolarization error randomly depolarizes in all basis angles at approximately the same rate, meaning consistent error in all cases. Note here that the overall error introduced by the incoherent noise model is much larger than the error introduced by the coherent noise model.  \n",
    "\n",
    "The errorless case, as expected, has effectively zero error for all qubit pairs, The tiny error seen is caused by rounding errors during simulation.\n",
    "\n",
    "The fused case shows both the relatively higher error from the incoherent error case, and the higher variance from the coherent error case. Across all qubits, this case produces the highest error seen yet, which makes sense given that both error models are in use simultaneously. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "c47ti1_uKnLs"
   },
   "source": [
    "### Circuit Fidelity by Cycle\n",
    "The qubit pair heatmaps are good, but summarize the fidelity data too much to reveal meaningful patterns. To rectify this, graph the circuit fidelity by cycle, alongside an exponential decay curve from the by-cycle fidelity estimation. For reference, an exponential decay associated with the original depolarization chance of `5e-3` is included. \n",
    "\n",
    "First some graphing utilities:\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "sL1nhnbyL_w9"
   },
   "outputs": [],
   "source": [
    "import seaborn as sns\n",
    "\n",
    "# map colors to qubit pairs\n",
    "colors = sns.cubehelix_palette(n_colors=graph.number_of_edges())\n",
    "colors = dict(zip(graph.edges, colors))\n",
    "\n",
    "\n",
    "def _p_gen(ax, rescale=False, decay=5e-3, color=None, colors=colors, label=None):\n",
    "    \"\"\"Create a _p function to graph circuit fidelities by layer.\n",
    "    It encodes the visualization options:\n",
    "    - rescale: rescale data under the exponential decay line to range [0,1].\n",
    "    - decay: exponent defining the exponential decay reference line.\n",
    "    - color: plot everything in given color, use by-qubit color map otherwise.\n",
    "    - colors: the by-qubit color map.\n",
    "    - label: label to give the data in the legend.\n",
    "    \"\"\"\n",
    "\n",
    "    def _p(record):\n",
    "        \"\"\"For a row of a fidelity data frame, plot based on encoded options\"\"\"\n",
    "\n",
    "        # pull fidelities from data frame row/record\n",
    "        data = record.fidelities\n",
    "\n",
    "        # plot horizontal at 1.0 fidelity\n",
    "        ax.axhline(1, color='grey', ls='--')\n",
    "\n",
    "        # compute exponential decay line from record parameters (a, layer_fid)\n",
    "        # these were computed by previous fitting from fidelity data\n",
    "        xx = np.linspace(0, record.cycle_depths.max())\n",
    "        exp = xeb_fitting.exponential_decay(xx, a=record.a, layer_fid=record.layer_fid)\n",
    "\n",
    "        # rescale data to expand area under exponential reference decay line\n",
    "        # the range from zero to the reference line becomes the range [0, 1]\n",
    "        # also plot that reference line in blue, which, once rescaled, is horizontal\n",
    "        if rescale:\n",
    "            data = data / ((1 - decay) ** (4 * record.cycle_depths))\n",
    "            exp = exp / ((1 - decay) ** (4 * xx))\n",
    "            ax.plot(xx, xx**0, label=r'Normalized Reference', color='blue')\n",
    "        else:\n",
    "            ax.plot(xx, (1 - decay) ** (4 * xx), label=r'Exponential Reference', color='blue')\n",
    "\n",
    "        _, _, pair = record.name\n",
    "        q0, q1 = pair\n",
    "        # use qubit pair to color map if color is not supplied, else just use color\n",
    "        plot_color = colors[pair] if color is None else color\n",
    "        # use qubit pair as label if not supplied, else just use label\n",
    "        plot_label = f'{q0}-{q1}' if label is None else label\n",
    "\n",
    "        # plot the data with the selected color and label\n",
    "        ax.plot(record.cycle_depths, data, 'o-', color=plot_color, label=plot_label)\n",
    "        # also plot the data-fitted exponential line\n",
    "        ax.plot(xx, exp, color=plot_color, alpha=0.3)\n",
    "\n",
    "    return _p\n",
    "\n",
    "\n",
    "def extras(ax, title):\n",
    "    \"\"\"Add labels, title and legend to an axis\"\"\"\n",
    "\n",
    "    ax.set_ylabel('Circuit fidelity')\n",
    "    ax.set_xlabel('Cycle Depth $d$')\n",
    "    ax.set_title(title)\n",
    "    # create a legend while making sure the handles are unique\n",
    "    handles, labels = ax.get_legend_handles_labels()\n",
    "    legend_items = dict(zip(labels, handles))\n",
    "    ax.legend(legend_items.values(), legend_items.keys(), loc='best')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "GlQSeoQFMFBy"
   },
   "source": [
    "Then use them to graph the different datasets."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "8c08c9ab8109"
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 2, figsize=(18, 13))\n",
    "for fidelities, ax, title in zip(fitted_fidelities, axes.flat, noise_titles):\n",
    "    fidelities.apply(_p_gen(ax), axis=1)\n",
    "    extras(ax, title)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "RcwjKxa7gGpB"
   },
   "source": [
    "The first, coherent error case, once again demonstrates the significant variance in circuit fidelity across the different qubits. This time it also shows that the variance between qubits can itself vary with the number of cycles that have been executed, which wasn't visible in the per-cycle estimation from the heatmaps. Additionally, individual exponential decay curves are fit and shown faded out for each qubit. The data only loosely fits an exponential decay, so using an exponential decay to model this type of coherent gate perturbation error may be imprecise. \n",
    "\n",
    "For the incoherent depolarization noise data, very little variance is seen between qubits, and all of them very closely match the reference curve. This empirically reinforces the expectation that a pauli depolarization noise model produces an exponential decay in the overall circuit fidelity. \n",
    "\n",
    "Unsurprisingly, the noiseless data produces no decrease in fidelity. \n",
    "\n",
    "In the fused case, both coherent perturbation and incoherent depolarization error contribute to the lowest circuit fidelity out of all cases. The overall fidelity is noticeably worse than the exponential reference and the incoherent data, as if they've been added together. It shows meaningfully more variance between qubits than the incoherent data, and this variation somewhat follows the pattern seen in the coherent data, with more variance around cycles 40 and 80, and less near cycles 20 and 60. However, it's somewhat hard to see this pattern given that the data is squashed into the space below the exponential reference curve.\n",
    "\n",
    "In the following plot, the fused data is shown, rescaled under the expected incoherent error's `5e-3` exponential decay, alongside the original coherent data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "2APLgj-XWkVF"
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(18, 6.5), sharey=True)\n",
    "two_fidelities = [fitted_fidelities[0], fitted_fidelities[3]]\n",
    "titles = [r'Unscaled Coherent', r'Rescaled Fused']\n",
    "for fidelities, rescale, ax, title in zip(two_fidelities, [False, True], axes.flat, titles):\n",
    "    fidelities.apply(_p_gen(ax, rescale=rescale), axis=1)\n",
    "    extras(ax, title)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "Yig2-glhWlOp"
   },
   "source": [
    "The two graphs are almost identical, with the rescaled fused data demonstrating the same patterns of by-qubit and by-cycle variation as the coherent-only data. It seems that the rescaling has approximately eliminated the influence of the incoherent depolarization error portion from the fused error, leaving only the coherent perturbation error. \n",
    "\n",
    "However, the two graphs are not exactly identical; there are small differences in the measured fidelities. Remember, both of these datasets were generated with identical circuits, random seeds, and other parameters. The noiseless fidelities plot from before demonstrated that rounding errors produce extremely small variances in fidelities. Even scaled up, it seems as though rounding errors are not sufficient to explain the difference between the coherent and rescaled-fused data. Instead, it is likely that the incoherent and coherent error included in the fused fidelity data are not entirely independent of one another in these benchmarking circuits. They likely interact nontrivially, producing some overlapping or unique error cases that would be found by, respectively, both or neither error model in isolation. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6025a292d19b"
   },
   "source": [
    "## Characterize Gate Parameters\n",
    "\n",
    "A classical optimizer can be used to approximate the two-qubit gate parameters that were actually used from the fidelity data. Here, the parameter of interest is the `phi` parameter, which is the one used to perturb `SQRT_ISWAP`. To characterize `phi`, flag it in the `SqrtISwapXEBOptions` object. Don't characterize the others parameters, for runtime reasons. This object, once passed into the optimizer wrapper `characterize_phased_fsim_parameters_with_xeb_by_pair`, denotes that the optimizer expects the two-qubit entangling gate used to be a `SQRT_ISWAP`, and that it will look for parameters to that gate (or the PhasedFSim version of it) which most closely replicate the observed fidelity data. \n",
    "\n",
    "In a real experiment, there are certainly additional kinds of coherent error. The following optimization procedure attempts to identify coherent error that produces a consistent change in the parameters of the circuits' two-qubit entangling gate, assuming it is in the `PhasedFSim` family. Other types of error, coherent and otherwise, may confound this process. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "baff45b4ad70"
   },
   "outputs": [],
   "source": [
    "import multiprocessing\n",
    "\n",
    "pool = multiprocessing.get_context('spawn').Pool()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "a3aaf1dfce4e"
   },
   "outputs": [],
   "source": [
    "# Set which angles to characterize (all)\n",
    "options = xeb_fitting.SqrtISwapXEBOptions(\n",
    "    characterize_theta=False,\n",
    "    characterize_zeta=False,\n",
    "    characterize_chi=False,\n",
    "    characterize_gamma=False,\n",
    "    characterize_phi=True,\n",
    ")\n",
    "# Parameterize the sqrt(iswap)s in circuit library\n",
    "pcircuits = [xeb_fitting.parameterize_circuit(circuit, options) for circuit in circuit_library]\n",
    "\n",
    "\n",
    "def characterization_result(sampled_df):\n",
    "    # Run the characterization loop\n",
    "    return xeb_fitting.characterize_phased_fsim_parameters_with_xeb_by_pair(\n",
    "        sampled_df,\n",
    "        pcircuits,\n",
    "        cycle_depths,\n",
    "        options,\n",
    "        pool=pool,\n",
    "        # ease tolerance so it converges faster:\n",
    "        fatol=1e-2,\n",
    "        xatol=1e-2,\n",
    "    )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "o-3iEEmIdcfl"
   },
   "source": [
    "Perform the coherent error characterization on the sample data, across all used noise models.\n",
    "\n",
    "Note: This cell will take 8 or so minutes to run on a default Google Colab instance. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "6-N3xeb8lN52"
   },
   "outputs": [],
   "source": [
    "%%time\n",
    "characterization_results = list(map(characterization_result, sampled_dfs))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "lc5QV4a8dkTI"
   },
   "source": [
    "Inspect the characterized parameters estimated for the coherent-only perturbation noise model. Note that, across all qubits, the characterization for the first, coherent error case, has correctly estimated the `phi` parameter to be in the range $[0, 2/16*\\pi] = [0, 0.393]$, or very close to it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "hZIvmHmMleoB"
   },
   "outputs": [],
   "source": [
    "characterization_results[0].final_params"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "6491ece85617"
   },
   "source": [
    "By additionally inspecting the original `phi_angles` used to add noise, you can see that the characterization was able to correctly estimate each qubit's individual phi value within about $\\pm 0.04$ error."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "de44f9b05329"
   },
   "outputs": [],
   "source": [
    "for qubits, phi in phi_angles.items():\n",
    "    print(qubits, phi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "qvXyscaYfpUv"
   },
   "source": [
    "## Visualizing Refitted Fidelities by Noise Model\n",
    "\n",
    "Finally, plot the original and refit circuit fidelities by cycle, with exponential decay curves, for each of the available datasets. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "1_Jfq4P5Vfll"
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(2, 2, figsize=(18, 13))\n",
    "result_dfs = list(map(lambda x: x.fidelities_df, characterization_results))\n",
    "result_fidelities = list(map(xeb_fitting.fit_exponential_decays, result_dfs))\n",
    "before_dfs = fitted_fidelities\n",
    "after_dfs = result_fidelities\n",
    "for before_df, after_df, ax, title in zip(before_dfs, after_dfs, axes.flat, noise_titles):\n",
    "    before_df.apply(_p_gen(ax, color='green', label=r'Original Fidelities'), axis=1)\n",
    "    after_df.apply(_p_gen(ax, color='orange', label=r'Refit Fidelities'), axis=1)\n",
    "    extras(ax, title)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "E4xIRDi9gBnu"
   },
   "source": [
    "In the coherent error case, the optimization was able to refit the fidelities very well. The new curves of the refit data almost perfectly match the horizontal line at `1.0` fidelity. The success of this refit means that the estimated parameters for this case are very likely to represent the actual parameters of the gate operation used, which was seen in the fact that the optimization re-discovered the `phi` value that was originally used to perturb `SQRT_ISWAP`. In the real hardware case, identifying these parameters means identifying what unitary the gate control hardware actually implemented. If this unitary is far from the intended one, hardware re-calibration or compensation must be performed in order to improve circuit fidelity. \n",
    "\n",
    "In constrast, the optimizer was completely unable to refit the circuit fidelities in the incoherent error case. The refit fidelities find no meaningful fidelity improvement, instead following the same, expected `5e-3` exponential decay as the original fidelities. Unsurprisingly, the optimizer was unable to find any coherent error in the incoherent-only fidelity data. \n",
    "\n",
    "The noiseless case also fails, but in a different way. Due to the fact that all of the original fidelities are extremely close to `1.0`, the optimization likely overvalues the miniscule differences in rounding error and overfits its model. \n",
    "\n",
    "The fused error case produces refit fidelities that are noticeably improved over the original ones, but only reach up to the exponential reference curve. It seems as though the optimization was able to detect and refit for the coherent error, but not the incoherent error. This means this parameter fitting procedure is still able to find improvements when multiple types of error are present. However, the fused graph again compacts all of the data under the exponential reference curve. The next cell performs the same rescaling as before to inspect the coherent error present in the fused data. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "id": "6Y4uDVpBry6B"
   },
   "outputs": [],
   "source": [
    "fig, axes = plt.subplots(1, 2, figsize=(18, 6.5), sharey=True)\n",
    "before_dfs = [fitted_fidelities[0], fitted_fidelities[3]]\n",
    "after_dfs = [result_fidelities[0], result_fidelities[3]]\n",
    "titles = [r'Unscaled Coherent', r'Rescaled Fused']\n",
    "for before_df, after_df, ax, title, rescale in zip(\n",
    "    before_dfs, after_dfs, axes.flat, titles, [False, True]\n",
    "):\n",
    "    before_df.apply(\n",
    "        _p_gen(ax, color='green', label=r'Original Fidelities', rescale=rescale), axis=1\n",
    "    )\n",
    "    after_df.apply(_p_gen(ax, color='orange', label=r'Refit Fidelities', rescale=rescale), axis=1)\n",
    "    extras(ax, title)\n",
    "\n",
    "fig.tight_layout()\n",
    "fig.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "utARIXLgrzNh"
   },
   "source": [
    "As noted, the optimizer is able to nearly perfectly refit the fidelities in the unscaled coherent error case. \n",
    "\n",
    "Rescaling the fused error data provides a similar result in the refit fidelities as in the original ones. The optimizer was unable to produce refit data with curves that follow the constant `1.0` line (which is the `5e-3` exponential decay curve when unscaled). Instead, the variety of orange refit curves around the blue reference line indicate that the presence of the incoherent error acts as a confounding variable to the optimizer, reducing its effectiveness noticeably. \n",
    "\n",
    "However, it still performed quite well overall, finding refit fidelities that are significantly better than the original ones. Additionally, the patterns of variance between qubits are similar between the refit fidelities of the coherent and rescaled data. This may imply that the optimizer is finding similar patterns of coherent error in the fused data as in the coherent-only data. Finally, the optimizer was able to identify the coherent error, even though the incoherent error was roughly an order of magnitude larger. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "E8BdVyNGLGnm"
   },
   "source": [
    "## Conclusion\n",
    "\n",
    "Cross-Entropy Benchmarking and optimizer refitting has been shown here to effectively characterize patterns of coherent error, to find the (parameters of the) true unitary operation used on individual qubit pairs. Importantly, this is effective even in the case where incoherent error is also acting on the system, but with noticeably reduced accuracy. In a real hardware system, with many interacting sources of error, XEB can still be useful to identify consistent coherent error, but it's important to remember that other error can confound these results. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "id": "wxvZ2plRMB-l"
   },
   "source": [
    "## What's Next?\n",
    "\n",
    "Now that you've identified coherent error, what can you do about it? \n",
    "- Depending on the parameter(s) identified to be out of tune, you may be able to change the circuits you want to run to compensate for that change.\n",
    "- In the real hardware case, the device may need to be re-calibrated in order to fix the tuning of particular gate operations on particular qubits. Reach out to your Google contact to let us know!"
   ]
  }
 ],
 "metadata": {
  "colab": {
   "collapsed_sections": [],
   "name": "coherent_vs_incoherent_xeb.ipynb",
   "toc_visible": true
  },
  "kernelspec": {
   "display_name": "Python 3",
   "name": "python3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
